# Replication files for:
# B�ck, Hanna, Marc Debus, and Michael Imre. 2022. "Populist Radical Parties, Pariahs, and Coalition Bargaining Delays." Party Politics.

### Part 1: Setup --------------------------------------------------------------
# R version used: 4.2.1
rm(list=ls())
options(scipen=10, digits=4)
mywd <- ""
setwd(mywd)

# install packages (if not installed yet) and load them
p_needed <- c("VGAM","MASS","survival","corpcor",
              "ggplot2","stargazer","dataverse")
# package versions used: corpcor_1.6.10, VGAM_1.1-5, dataverse_0.3.10, stargazer_5.2.2, survival_3.2-13, ggplot2_3.3.6, MASS_7.3-55     
packages <- rownames(installed.packages())
p_to_install <- p_needed[!(p_needed %in% packages)]
if (length(p_to_install) > 0) {
  install.packages(p_to_install)
}

lapply(p_needed, require, character.only = TRUE)

rm(list=ls())

# load data
load("baeck_et_al_2022_data.RData")

## Variables included in the data (more thorough definitions for the variables used in the models can be found in the article)
# Bundesland: state
# year: year
# election: election
# choice: binary; was there a bargaining attempt with these parties?
# failed: binary; did the bargaining attempt fail?
# govform_duration: number of days the bargaining attempt took
# cdu-afd: binary variables indicating whether the specific party is part of the respective potential attempt
# inc: Incumbent coalition
# ccoal: Cross-cutting coalition
# strongest: Largest party in parliament included
# phet: Intra-cabinet policy heterogeneity
# rext_s: Seat share right-wing extremist parties
# lext_s: Seat share left-wing extremist parties
# pariah: Pariah party among bargaining parties
# lr_esteban_ray: Ideological polarization in parliament
# blocmajority: Majority for centre-left or centre-right bloc
# enp: Effective number of parliamentary parties
# single_partygov: Single party majority government
# decadeXXXX: decade dummies
# Bundesland_XXXX: state dummies


### Part 2: Descriptive figures ------------------------------------------------
## Figure 1
tiff("figure1.tif", res=600, compression = "lzw", height=7, width=10, units="in")
hist(data$govform_duration, breaks = seq(0,300,20), border=T,
     main="",xlab="Duration (in days)",ylim=c(0,50),las=1,col="grey70")
dev.off()


## Figure 2
fig2data <- data[data$choice == 1 | data$failed==1,]

tiff("figure2.tif", res=600, compression = "lzw", height=7, width=10, units="in")
ggplot(fig2data, aes(x=as.factor(Bundesland), y=lr_esteban_ray)) + 
  geom_boxplot(fill="grey70") + 
  xlab("") +  ylab("Ideological polarisation within parliaments") + theme_bw() +
  theme(axis.text = element_text(size = 11.5)) +
  scale_y_continuous(breaks=seq(0,3,0.5)) +
  theme(axis.text.x=element_text(angle=55,hjust=1))
dev.off()

### Part 3: Main Models --------------------------------------------------------
# this whole section, to a large degree, closely follows the approach outlined in the replication materials of Ecker and Meyer (2020)

# load log-likelihood functions defined by Chiba, Martin, and Stevenson (2015)
# adjusted for updated version of VNORM function by Ecker and Meyer (2020)
# check if the script file from the replication materials of Ecker and Meyer (2020) already exists in Folder, if not download it from Dataverse
if (!file.exists("functions_ecker_meyer_2020.R")) {
  writeBin(get_file("3_LLfunctions_newVNORM.R", 
                    "doi:10.7910/DVN/2XSNPC", 
                    server = "dataverse.harvard.edu"),"functions_ecker_meyer_2020.R")
}

source("functions_ecker_meyer_2020.R")


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
### Model part 1: formation stage

# create a variable called caseid (necessary as grouping variable in clogit model)
data$caseid <- data$election

# formation equation
eqn.formation <- as.formula(choice ~ inc + ccoal + strongest + phet)

# estimate formation stage
univ.formation <- estimate.clogit(eqn.formation, data, constraint = F)

# store estimation results
save(univ.formation, file = "univ.formation.rda")

# use estimates as starting values for joint estimation
beta.formation <- univ.formation$fit$par

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
### Model part 2: duration stage

### Model 2.1
# duration equation
eqn.duration1 <- as.formula(govform_duration ~ 
                              rext_s + lext_s +
                              lr_esteban_ray +
                              enp + single_partygov + 
                              decade2000 + decade2010 + decade2020 +
                              Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                              Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                              Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                              Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                              Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia)

# estimate duration stage
univ.duration1 <- estimate.univdur(eqn.duration1, "election", "weibull", data)

# store estimation results
save(univ.duration1, file = "univ.duration1.rda")

# use estimates as starting values for joint estimation
beta.duration1 <- univ.duration1$fit$par

# bivariate weibull analysis
# starting values
joint.init1 <- c(beta.formation, beta.duration1, 0)

# estimate selection-duration dependent competing-risk Weibull AFT model
biv.formation.duration1 <- estimate.cld(eqn.formation, eqn.duration1, 
                                       failure.type = "election", dist = "weibull", data = data,
                                       WantHessian = T, init = joint.init1, constraint = F)

# store estimation results
save(biv.formation.duration1, file = "biv.formation.duration1.rda")

# display empirical results Table 1
m1 <- disp.cld(biv.formation.duration1)
biv.formation.duration1$fit$value


### Model 2.2
# duration equation
eqn.duration2 <- as.formula(govform_duration ~ 
                              pariah  +
                              lr_esteban_ray +
                              enp + single_partygov + 
                              decade2000 + decade2010 + decade2020 +
                              Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                              Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                              Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                              Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                              Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia)

# estimate duration stage
univ.duration2 <- estimate.univdur(eqn.duration2, "election", "weibull", data)

# store estimation results
save(univ.duration2, file = "univ.duration2.rda")

# use estimates as starting values for joint estimation
beta.duration2 <- univ.duration2$fit$par

# bivariate weibull analysis
# starting values
joint.init2 <- c(beta.formation, beta.duration2, 0)

# estimate selection-duration dependent competing-risk Weibull AFT model
biv.formation.duration2 <- estimate.cld(eqn.formation, eqn.duration2, 
                                        failure.type = "election", dist = "weibull", data = data,
                                        WantHessian = T, init = joint.init2, constraint = F)

# store estimation results
save(biv.formation.duration2, file = "biv.formation.duration2.rda")

# display empirical results Table 1
m2 <- disp.cld(biv.formation.duration2)
biv.formation.duration2$fit$value


### Model 2.3
# duration equation
eqn.duration3 <- as.formula(govform_duration ~ 
                              rext_s + lext_s + pariah  +
                              lr_esteban_ray +
                              enp + single_partygov + 
                              decade2000 + decade2010 + decade2020 + 
                              Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                              Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                              Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                              Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                              Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia)

# estimate duration stage
univ.duration3 <- estimate.univdur(eqn.duration3, "election", "weibull", data)

# store estimation results
save(univ.duration3, file = "univ.duration3.rda")

# use estimates as starting values for joint estimation
beta.duration3 <- univ.duration3$fit$par

# bivariate weibull analysis
# starting values
joint.init3 <- c(beta.formation, beta.duration3, 0)

# estimate selection-duration dependent competing-risk Weibull AFT model
biv.formation.duration3 <- estimate.cld(eqn.formation, eqn.duration3, 
                                        failure.type = "election", dist = "weibull", data = data,
                                        WantHessian = T, init = joint.init3, constraint = F)

# store estimation results
save(biv.formation.duration3, file = "biv.formation.duration3.rda")

# display empirical results Table 1
m3 <- disp.cld(biv.formation.duration3)
biv.formation.duration3$fit$value

# export all models
stargazer(m1,m2,m3,type="text",out="table1.doc")



### Part 4: Plotting effect sizes ----------------------------------------------
# code in this section was mostly adapted from the replication materials of Chiba et al. (2015)

## subset the data
dur.pool <- data[data$choice==1,]
dur.pool.r <- dur.pool[dur.pool$failed==0,]

## Range of observed duration until replacement
range(dur.pool.r$govform_duration)
summary(dur.pool.r$govform_duration)
summary(dur.pool$govform_duration)

ft <- seq(from = 1, to = 281, length = 281)

## make an empty dataframe for all setx
setxdataframe <- data.frame(matrix(NA, nrow = 14, ncol = 25))
## make an empty dataframe for all linedata
linesdataframe <- data.frame(matrix(NA, nrow = 14, ncol = 281))

# set data for all scenarios
setxdataframe[1,] <- c(median(dur.pool$rext_s), median(dur.pool$lext_s),
              1, mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[2,] <- c(mean(dur.pool$rext_s), mean(dur.pool$lext_s),
              0, mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[3,] <- c(0, mean(dur.pool$lext_s), 
              median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[4,] <- c(0.05, mean(dur.pool$lext_s), 
              median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[5,] <- c(0.1, mean(dur.pool$lext_s), 
              median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[6,] <- c(0.15, mean(dur.pool$lext_s), 
              median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[7,] <- c(0.2, mean(dur.pool$lext_s),
              median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[8,] <- c(0.25, mean(dur.pool$lext_s), 
              median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[9,] <- c(mean(dur.pool$rext_s), 0,
              median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[10,] <- c(mean(dur.pool$rext_s), 0.05, 
              median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[11,] <- c(mean(dur.pool$rext_s), 0.1,
              median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
              mean(dur.pool$enp), median(dur.pool$single_partygov),
              rep(0,18),  # fixed effects
              1) # constant

setxdataframe[12,] <- c(mean(dur.pool$rext_s), 0.15, 
               median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
               mean(dur.pool$enp), median(dur.pool$single_partygov),
               rep(0,18),  # fixed effects
               1) # constant

setxdataframe[13,] <- c(mean(dur.pool$rext_s), 0.2, 
               median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
               mean(dur.pool$enp), median(dur.pool$single_partygov),
               rep(0,18),  # fixed effects
               1) # constant

setxdataframe[14,] <- c(mean(dur.pool$rext_s), 0.25, 
               median(dur.pool$pariah), mean(dur.pool$lr_esteban_ray),
               mean(dur.pool$enp), median(dur.pool$single_partygov),
               rep(0,18),  # fixed effects
               1) # constant

# create estimates for all scenarios
for (x in 1:14) {
  setx <- as.numeric(setxdataframe[x,])
  
  ## Survivor function from Joint Model
  nk <- length(biv.formation.duration3$fit$par)
  cw.p <- exp(biv.formation.duration3$fit$par[(nk-1)])
  cw.b <- biv.formation.duration3$fit$par[5:29]
  cw.l <- exp( -setx %*% cw.b )
  cw.St <- exp(-(ft*cw.l)^cw.p)
  cw.ft <- (cw.p*cw.l) * (ft*cw.l)^(cw.p-1) * exp(- (ft*cw.l)^cw.p)
  
  ## Median duration
  cw.l^(-1) * log(2)^(1/cw.p)
  
  setx.L <- c(
    quantile(dur.pool$rext_s)[2],
    quantile(dur.pool$lext_s)[2],
    quantile(dur.pool$pariah)[2], mean(dur.pool$lr_esteban_ray),
    mean(dur.pool$enp), median(dur.pool$single_partygov), rep(0,18), 1)
  
  setx.H <- c(
    quantile(dur.pool$rext_s)[4],
    quantile(dur.pool$lext_s)[4],
    quantile(dur.pool$pariah)[4], mean(dur.pool$lr_esteban_ray),
    mean(dur.pool$enp), median(dur.pool$single_partygov), rep(0,18), 1)
  
  ## FD in survivor function from Joint Model
  cw.l.L <- exp(-setx.L %*% cw.b)
  cw.St.L <- exp(-(ft * cw.l.L)^cw.p)
  cw.l.H <- exp(-setx.H %*% cw.b)
  cw.St.H <- exp(-(ft * cw.l.H)^cw.p)
  cw.fd <- cw.St.H - cw.St.L
  
  ## First difference in terms of median duration
  cw.l.H^(-1) * log(2)^(1/cw.p) - cw.l.L^(-1) * log(2)^(1/cw.p)
  
  # uncertainty estimates
  n.sim <- 1000
  nw.fd.ci <- cw.fd.ci <- nw.St.ci <- cw.St.ci <- matrix(NA, nrow=n.sim, ncol=length(ft))
  
  ## Estimated beta and v-cov matrices
  cw.beta <- biv.formation.duration3$fit$par
  cw.vcov <- -solve(biv.formation.duration3$fit$hessian)
  
  ## simulate Beta
  set.seed(123456)
  cw.simb <- mvrnorm(n.sim, cw.beta, cw.vcov)
  
  i <- 1
  for (i in 1:n.sim){
    ## Betas for this round of simulation
    cw.p <- exp(cw.simb[i, (nk-1)])
    cw.b <- cw.simb[i,c(5:29)]
    
    ## Survivor function from Joint Model
    cw.l <- exp( -setx %*% cw.b)
    cw.St.ci[i,] <- exp(-(ft * cw.l)^cw.p)
    
    ## FD in survivor function from Joint Model
    cw.l.L <- exp( -setx.L %*% cw.b)
    cw.St.L <- exp(-(ft * cw.l.L)^cw.p)
    cw.l.H <- exp( -setx.H %*% cw.b)
    cw.St.H <- exp(-(ft * cw.l.H)^cw.p)
    cw.fd <- cw.St.H - cw.St.L
    cw.fd.ci[i,] <- cw.fd
  }
  
  cw.St.q <- apply(cw.St.ci, MARGIN = 2, FUN = quantile, probs = c(.025, .5, .975))
  cw.fd.q <- apply(cw.fd.ci, MARGIN = 2, FUN = quantile, probs = c(.025, .5, .975))

  linesdataframe[x,] <- cw.St.q[2,]
}
  
x.rng <- c(0, 285); y.rng <- c(0, 1)
par(mar=c(2.5,3,2,2)+0.1, mgp=c(1.5,0.5,0), bg="white",
    cex = 1, cex.main = 1, cex.lab=1,cex.axis=1)

## Figure 3
tiff("figure3.tif", res=600, compression = "lzw", height=7, width=10, units="in")
plot(x = x.rng, y = y.rng, xaxs="i", yaxs="i", main = "", type = "n",
     xlab = "Time (in days)", ylab = "Survival Probability",axes=F)
box(); axis(2)
axis(1, at=seq(from=0,to=300, by=20), labels=seq(from=0,to=300, by=20))
abline(h=0, col="gray80")

lines(ft, linesdataframe[3,], col="black", lwd=2, lty=1)
lines(ft, linesdataframe[4,], col="gray20", lwd=2, lty=5)
lines(ft, linesdataframe[5,], col="gray40", lwd=2, lty=6)
lines(ft, linesdataframe[6,], col="gray60", lwd=2, lty=4)
lines(ft, linesdataframe[7,], col="black", lwd=2, lty=2)
lines(ft, linesdataframe[8,], col="gray20", lwd=2, lty=3)

legend(x = "topright", lty = c(1,5,6,4,2,3), text.font = 1, 
       col= c("black","gray20","gray40","gray60","black","gray20"),text.col = "black", 
       legend=c("Seat share RRP=0.00", "Seat share RRP=0.05", "Seat share RRP=0.10", 
                "Seat share RRP=0.15", "Seat share RRP=0.20", "Seat share RRP=0.25"))

dev.off() 


## Figure 4
tiff("figure4.tif", res=600, compression = "lzw", height=7, width=10, units="in")
plot(x = x.rng, y = y.rng, xaxs="i", yaxs="i", main = "", type = "n",
     xlab = "Time (in days)", ylab = "Survival Probability",axes=F)
box(); axis(2)
axis(1, at=seq(from=0,to=300, by=20), labels=seq(from=0,to=300, by=20))
abline(h=0, col="gray80")

lines(ft, linesdataframe[9,], col="black", lwd=2, lty=1)
lines(ft, linesdataframe[10,], col="gray20", lwd=2, lty=5)
lines(ft, linesdataframe[11,], col="gray40", lwd=2, lty=6)
lines(ft, linesdataframe[12,], col="gray60", lwd=2, lty=4)
lines(ft, linesdataframe[13,], col="black", lwd=2, lty=2)
lines(ft, linesdataframe[14,], col="gray20", lwd=2, lty=3)

legend(x = "topright", lty = c(1,5,6,4,2,3), text.font = 1, 
       col= c("black","gray20","gray40","gray60","black","gray20"),text.col = "black", 
       legend=c("Seat share RLP=0.00", "Seat share RLP=0.05", "Seat share RLP=0.10", 
                "Seat share RLP=0.15", "Seat share RLP=0.20", "Seat share RLP=0.25"))

dev.off()


## Figure 5
tiff("figure5.tif", res=600, compression = "lzw", height=7, width=10, units="in")
plot(x = x.rng, y = y.rng, xaxs="i", yaxs="i", main = "", type = "n",
     xlab = "Time (in days)", ylab = "Survival Probability",axes=F)
box(); axis(2)
axis(1, at=seq(from=0,to=300, by=20), labels=seq(from=0,to=300, by=20))
abline(h=0, col="gray80")

lines(ft, linesdataframe[1,], col="black", lwd=2, lty=1)
lines(ft, linesdataframe[2,], col="gray50", lwd=2, lty=2)

legend(x = "topright", lty = c(1,2), text.font = 1, 
       col= c("black","gray50"),text.col = "black", 
       legend=c("Pariah party", "No pariah party"))
dev.off() 



### Part 5: Alternative model specifications (Appendices) ----------------------

# Appendix A ----
# duration equation
eqn.duration4 <- as.formula(govform_duration ~ 
                              rext_s + lext_s + pariah  +
                              lr_esteban_ray + blocmajority + blocXestebanray +
                              enp + single_partygov +
                              decade2000 + decade2010 + decade2020 + 
                              Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                              Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                              Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                              Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                              Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia)

# estimate duration stage
univ.duration4 <- estimate.univdur(eqn.duration4, "election", "weibull", data)

# store estimation results
save(univ.duration4, file = "univ.duration4.rda")

# use estimates as starting values for joint estimation
beta.duration4 <- univ.duration4$fit$par

# bivariate weibull analysis
# starting values
joint.init4 <- c(beta.formation, beta.duration4, 0)

# estimate selection-duration dependent competing-risk Weibull AFT model
biv.formation.duration4 <- estimate.cld(eqn.formation, eqn.duration4, 
                                        failure.type = "election", dist = "weibull", data = data,
                                        WantHessian = T, init = joint.init4, constraint = F)

# store estimation results
save(biv.formation.duration4, file = "biv.formation.duration4.rda")

# display empirical results Table 1
m4 <- disp.cld(biv.formation.duration4)
biv.formation.duration4$fit$value

# export all models
stargazer(m4,type="text",out="tableA1.doc")


# Appendix B ----
model1.weibull <- survreg(Surv(govform_duration,censor_discelec)~ 
                              rext_s + lext_s +
                              lr_esteban_ray +
                              enp + single_partygov + 
                              decade2000 + decade2010 + decade2020 +
                              Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                              Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                              Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                              Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                              Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia, data=data,dist="weibull")

model2.weibull <- survreg(Surv(govform_duration,censor_discelec)~ 
                            pariah  +
                            lr_esteban_ray +
                            enp + single_partygov + 
                            decade2000 + decade2010 + decade2020 +
                            Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                            Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                            Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                            Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                            Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia, data=data,dist="weibull")

model3.weibull <- survreg(Surv(govform_duration,censor_discelec)~ 
                            rext_s + lext_s + pariah  +
                            lr_esteban_ray +
                            enp + single_partygov + 
                            decade2000 + decade2010 + decade2020 + 
                            Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                            Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                            Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                            Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                            Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia, data=data,dist="weibull")

model4.weibull <- survreg(Surv(govform_duration,censor_discelec)~ 
                            rext_s + lext_s + pariah  +
                            lr_esteban_ray + blocmajority + blocXestebanray +
                            enp + single_partygov +
                            decade2000 + decade2010 + decade2020 + 
                            Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                            Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                            Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                            Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                            Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia, data=data,dist="weibull")

stargazer(model1.weibull,model2.weibull,model3.weibull,model4.weibull,
          type="html",out="tableA2.html",star.char = c("+", "*", "**"),
          star.cutoffs = c(0.1, 0.05, 0.01), 
          column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),
          model.numbers=FALSE, dep.var.caption  = "",dep.var.labels.include=FALSE,
          keep=c("rext_s","lext_s","pariah","lr_esteban_ray","blocmajority",
                "blocXestebanray","enp","single_partygov"),
          covariate.labels=c("Seat share right-wing extremist parties","Seat share left-wing extremist parties",
                            "Pariah party among bargaining parties",
                             "Ideological polarization in parliament","Majority for centre-left or centre-right bloc",
                             "Maj. for bloc X Ideolog. polarization","Effective no. parliamentary parties",
                             "Single party government"),
          add.lines=list(c("State dummies included",rep("Yes",4)),
                         c("Decade dummies included",rep("Yes",4)),
                         c("Log(scale)", c(signif(log(model1.weibull$scale),3),
                                           signif(log(model3.weibull$scale),3),
                                           signif(log(model4.weibull$scale),3)))))


# Appendix C ----  
model1.cox <- coxph(Surv(govform_duration,censor_discelec)~ 
                        rext_s + lext_s +
                            lr_esteban_ray +
                            enp + single_partygov + 
                            decade2000 + decade2010 + decade2020 +
                            Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                            Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                            Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                            Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                            Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia, data=data)

model2.cox <- coxph(Surv(govform_duration,censor_discelec)~ 
                            pariah  +
                            lr_esteban_ray +
                            enp + single_partygov + 
                            decade2000 + decade2010 + decade2020 +
                            Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                            Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                            Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                            Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                            Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia, data=data)

model3.cox <- coxph(Surv(govform_duration,censor_discelec)~ 
                            rext_s + lext_s + pariah  +
                            lr_esteban_ray +
                            enp + single_partygov + 
                            decade2000 + decade2010 + decade2020 + 
                            Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                            Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                            Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                            Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                            Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia, data=data)

model4.cox <- coxph(Surv(govform_duration,censor_discelec)~ 
                            rext_s + lext_s + pariah  +
                            lr_esteban_ray + blocmajority + blocXestebanray +
                            enp + single_partygov +
                            decade2000 + decade2010 + decade2020 + 
                            Bundesland_Bavaria + Bundesland_Berlin + Bundesland_Brandenburg +      
                            Bundesland_Bremen + Bundesland_Hamburg + Bundesland_Hesse +        
                            Bundesland_Mecklenburg_West_Pommerania + Bundesland_Lower_Saxony + Bundesland_North_Rhine_Westphalia +        
                            Bundesland_Rhineland_Palatinate + Bundesland_Saarland + Bundesland_Saxony +            
                            Bundesland_Saxony_Anhalt + Bundesland_Schleswig_Holstein + Bundesland_Thuringia, data=data)


stargazer(model1.cox,model2.cox,model3.cox,model4.cox,
          type="html",out="tableA3.html",star.char = c("+", "*", "**"),
          star.cutoffs = c(0.1, 0.05, 0.01),          
          column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),
          model.numbers=FALSE, dep.var.caption  = "",dep.var.labels.include=FALSE,
          keep=c("rext_s","lext_s","pariah","lr_esteban_ray","blocmajority",
                 "blocXestebanray","enp","single_partygov"),
          covariate.labels=c("Seat share right-wing extremist parties","Seat share left-wing extremist parties",
                             "Pariah party among bargaining parties",
                             "Ideological polarization in parliament","Majority for centre-left or centre-right bloc",
                             "Maj. for bloc X Ideolog. polarization","Effective no. parliamentary parties",
                             "Single party government"),
          add.lines=list(c("State dummies included",rep("Yes",4)),
                         c("Decade dummies included",rep("Yes",4))))




### Part 6: Cleaning -----------------------------------------------------------

# delete files created during script but not needed as outputs
files <- c("univ.formation.RDA","univ.duration1.RDA","univ.duration2.RDA",
           "univ.duration3.RDA","univ.duration4.RDA","biv.formation.duration1.RDA",
           "biv.formation.duration2.RDA","biv.formation.duration3.RDA","biv.formation.duration4.RDA")
for (file in 1:length(files)) {
  if (file.exists(files[file])) {
    file.remove(files[file])
  }
}

